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Abstract 

We calculate the cosmic ray spectrum of ultra high energy neutralinos that 
one should expect provided that the observed ultra high energy cosmic rays 
are produced by the decay of superheavy particles X, Mx > 10 12 GeV, in 
supersymmetric models. Our calculation uses an extended DGLAP formal- 
ism. Forthcoming cosmic ray observatories should be able to detect these 
neutralinos. 
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I. INTRODUCTION 



Over one hundred cosmic ray events with energy higher than 4x 10 10 GeV have been detected by 
past and present observatories. About 30 events have energy in the range (1 — 3) x 10 11 GeV (see [|J 
for an experimental review). In top-down production mechanisms of Ultra High Energy Cosmic 
Rays (UHECR) the decay products of very massive particles X with Mx > 10 12 GeV account for 
these events. Particle X can be either a GUT scale particle produced by the decay of a cosmological 
network of topological defects created in a phase transition in the early universe 0-[§, or a cold 
dark matter particle clustered in the galactic halo whose abundance was generated during or shortly 
after inflation [|6| P~C|] - Whatever the nature of X is, it decays into partons which hadronise and give 
the UHECR primaries that we observe on the Earth. Therefore, in order to calculate the predicted 
spectrum and composition of UHECRs one has to calculate the nonperturbative Fragmentation 
Functions (FFs) of partons into hadrons at the energy scale Mx- Several approaches have been 
taken to solve the fragmentation process and calculate the spectrum of UHECRs: Montecarlo 
generators [§|[lj, analytical estimates (MLLA) @f|[l],[l2| and DGLAP evolution (see §Jl5] 

for a comparison of these different approaches). 

As long as supersymmetry (SUSY) is a low energy symmetry of nature, the decay of X will 
produce supersymmetric particles. If R-parity is conserved, the Lightest Supersymmetric Particle 
(LSP) is stable and should be among the primary cosmic rays that reach the Earth. Theoretical 
motivations favour the neutralino as the LSP, hence we shall concentrate on neutralinosQ. In the 
Minimal Supersymmetric Standard Model (MSSM) with R-parity conserved, neutralinos are a mix- 
ture of higgsinos and neutral gauginos (bino, B, and neutral wino, W°). The precise composition 
of the neutralinos in terms of the interaction eigenstates depends on the particular SUSY scenario. 
For simplicity, we shall concentrate on the Constrained MSSM, which assumes that the soft SUSY- 
breaking scalar masses, the gaugino masses and the trilinear parameters are all universal at the 
GUT scale. In this scenario, it can be shown that in the region of the parameter space where the 
relic density of neutralinos is of cosmological significance, the LSP is mainly a bino fl8[| . 

Weakly interacting particles like neutralinos or neutrinos have cross sections with ordinary 
matter too small to be detected with present day UHECR observatories. Baryons, nuclei and 



perhaps photons may account for all the events detected so far [19]. Forthcoming detectors will be 



sensitive to weakly interacting UHECRs (see Refs. [2^,21| and references therein). It is therefore 
important to obtain a good estimate of the neutralino flux expected in top-down models. 

The neutralino spectrum in top-down models has already been estimated using Montecarlo 



generators by Berezinsky&Kachelriefi [11,22|. In the Montecarlo approach a primary parton with 
energy Mx /2 produced in the decay of X initiates a parton cascade which proceeds until a specific 
minimal virtuality is reached. The squarks and gluinos in the parton shower decay into the LSP 
and ordinary partons once they reach the scale Msusy> the universal mass of squarks and gluinos. 
All quark and gluons evolve down to the hadronization scale where a phenomenological model is 
employed to bind them into hadrons. 

The Montecarlo simulations did not include the electroweak radiation of neutralinos by off- 
shell partons in the shower. In the next section we show how the DGLAP formalism employed to 



x In some SUSY models the gravitino, the gluino and even the sneutrino could be the LSP. An 
exotic hadron state containing a stable gluino has been proposed to explain the UHECR events [17]. 
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calculate the neutrino, photon and baryon spectra of UHECR |T^-^] can be extended to calculate 
the spectrum of neutralinos. We concentrate on the electroweak radiation of neutralinos by off- 
shell partons to complement the Montecarlo simulations. We find that at large x electroweak 
radiation of neutralinos is important. The large x region is the most relevant for future cosmic ray 
observatories since in this region the neutralino flux is larger than the neutrino flux and will not 
be swamped by the neutrino signal. 



II. EVOLUTION EQUATIONS FOR NEUTRALINO FRAGMENTATION 

FUNCTIONS 



The main channels of neutralino production are the hadronic decays of X. Hence the flux of 
neutralinos is given by the fragmentation functions from partons, D%(x, M x ), where a is any quark 
qk and squark flavour s^, a gluon g or a gluino A. The flux is proportional to the sum over a of 
D x (x, M x ) each one weighted with the branching ration for the decay of X into parton a. The 
variable < x < 1 is the fraction of the momentum of a carried off by \- These nonperturbative 
functions depend on the energy scale /u; in our case \x = Mx- The rate of change of D%(x, /j, 2 ) with 
/i 2 is given by the sum of two terms. The first one takes into account the ordinary SUSY QCD 
branching of partons and gives rise to the so-called Dokshitzer-Gribov-Lipatov-Altarelli-Parisi 
(DGLAP) equations [23,24|. The second term takes into account the emission of neutralinos in 
the parton shower, which stems from the tree level weak SUSY coupling among quarks, squarks 
and neutralinos. For any (anti)quark and (anti)squark we obtain, to lowest order in the strong 



coupling as and the weak couplings a^ k , 



^d^D^x^ 2 ) 
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P XSk (x)®DX(x,v 2 ), (2) 



where the dots inside the square brackets stand for the remaining leading order SUSY QCD terms^, 
see [[l5,25|,26|]. For gluons and gluinos, to lowest order in aw, there are no additional terms to the 
usual SUSY QCD ones. The FF of neutralinos from neutralinos is D*(x, fj?) = 5(1 — x) + 0(a\y)- 
Note that the functions D$(x, fi 2 ) are of order a^y . The splitting functions of quarks and squarks 



into neutralinos are given by 



\<ik 



1 — x and P, 



1. We make use of the convolution operator 



f 1 dz 

A{x)®B{x) = / — A(4B(§). 



(3) 



Analogous equations to Eqs. (|T(-|2t) employed to study photon structure or photon production by 



parton showers can be found in |27j2 



2 The DGLAP equations do not include soft gluon emission coherence. Coherence is important at 
low x. Our results will only hold for x much larger than the position of the Gaussian peak produced 
by coherence, at x p ~ J ' Kqqd/Mx- As pointed out in Sec. | we are only interested in the large x 
region. The DGLAP equations can be modified to include coherence, see for example [13]. 
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It is useful to introduce the following evolution variable 



2-Kb as{n 2 )'' 



(4) 



b being the coefficient in the leading order /3-function governing the running of the strong coupling: 
(3(a s ) = -ba 2 s . 

For simplicity we shall only consider gaugino production. We shall argue later that the final 
results are not substantially altered if higgsinos are included in the calculation. Also for clarity we 
shall concentrate on the evolution of quark and squark singlet functions, defined as the following 
sum over flavours 



k 



(5) 
(6) 



The singlet functions are coupled to the gluons and gluinos by means of the following 4x4 matrix 
integro-differential equation 
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(7) 



where = J2k anc ^ nF ^ s ^ ne num ber of flavours. These will be the only relevant equations 
if one assumes flavour universality in the decay of X. If one wishes to study the evolution of each 
flavour separately then one needs the complete set of DGLAP equations as given in HH^I to 
which one has to add the neutralino radiation terms on the right hand side of Eqs. 
We can write Eq. (|^) in the simple matrix form 



d r xD x = xP (8) xD x + xX x , 



(8) 



where xD x and xX x are 4— vectors and xP is a 4 x 4 matrix whose elements are all the leading 
order SUSY QCD splitting functions P a \, times x. These functions were calculated in [25,|FJ]. We 
have multiplied Eq. (|7]) by x to improve numerical stability. 

Equation @ is a linear inhomogeneous equation. Its associated homogeneous equation is the 
usual DGLAP equation for the coupled evolution of the quark and squark singlets, gluons and 
gluinos. The inhomogeneous or source term arises from the neutralino radiation by partons. The 
general solution to Eq. (§) is the sum of the general solution to the associated homogeneous equation 
plus a particular solution 



xD x (x, t) = E(x, r - t ) ® xD x (x, r ) + / dr'E{x, r - r') ® xX x (x, r). 



(9) 



TO 



initial condition E(x,0) = 5(1 — x)I± (see [16.29]) 



The evolution operator E(x, r) is the solution to the associated homogeneous equation with the 

(10) 



1 



E(x, t) = e xP{x)T = 5(1 - x)h + xP(x)t + ^xP(x) ® xP(x)t 2 + . 
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We take as initial energy scale hq = AfeuSY, to = T (^o) = 0, the typical scale for all sparticle 
masses. 

Equation (0) gives the total neutralino FF as the sum of two terms. The first one, the solution 
to the homogeneous equation, is the contribution from partons that decay on-shell into neutralinos 
(assuming a universal mass for all squarks and gluinos Mgusy)- This contribution was studied 



with Montecarlo generators in Refs. [11,22]. The second term is the contribution stemming from 
the electroweak radiation of neutralinos by partons. From now on, we concentrate on this second 
term and take 

xD x {x,t)= [ T dr'E(x,T -t')®xX x (x,t'). (11) 
Jo 

This expression is our main equation. It will allow us to estimate the relevance of electroweak 
radiative emission of neutralinos. 



III. NUMERICAL CALCULATION AND RESULTS 



In order to calculate numerically xD x (x,t) as given in Eq. ( |l l[) we expand the evolution 
operator in Laguerre polynomials^ L n (x) 

oo 

E(x,r) = J2 E n(r)L n (-\nx), (12) 

n=0 

4 n k 

£n(r)=£^ r £^>- (13) 

i=l k=0 

The scalars A« are the eigenvalues of xPq, the first coefficient in the Laguerre expansion of the 
matrix xP(x). Energy conservation gives Ai = 0. All other eigenvalues are negative. The 4x4 
matrices B\ n are given by recursive relations that were calculated in [16]. Likewise we expand the 
source term in Laguerre polynomials 

xX x (x,r) = ^lf(x), (14) 

CO 

f(x) = J2fnL n (-lnx), (15) 

n=0 



x ,27 3 \3J V2 

The Laguerre expansions for the neutralino FFs are 



;£>*(*, r) = £ C x (t)L n (-Iux), (17) 



N=0 



where Cq = Cq and = — C^_ 1 for N > 0. The 4— vectors Cjy are calculated substituting 
all Laguerre expansions in Eq. (p4|) 



3 Laguerre polynomial expansions to study scaling violations in QCD were first introduced in [29] 
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In order to present our results we multiply FFs by x 3 since the measured spectrum is usually 
multiplied by E 3 to highlight its structure. We plot in Fig. ([l]) our numerical calculation for the bino 
x 3 D^(x,M x ) when Mx = 10 12 GeV and Msusy = 400 GeV. The major contribution comes from 
squarks and quarks, which can radiate neutralinos at order aw- The contribution from gluinos 
and gluons is much smaller since they do not have order aw coupling to neutralinos, therefore 
to lowest order they can only generate neutralinos through mixing with squarks and quarks, see 
Eq. (0). The fraction of available energy Mx = 10 12 GeV carried by binos is (xD^(x,Mx) i 
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FIG. 1. B fragmentation functions for M x = 10 12 GeV and M SU SY = 400 GeV. The solid 
line is the squark contribution, the dotted lines is the quark one, the dashed line is the gluino 
contribution and the dot-dashed line is the gluon one. 



C^(t(M x )) = (0.047,0.002,0.109,0.003) for q,g, s,X, respectively (divide by 2n F to get fraction 
per quark and squark flavour). From Eqs. (|ll|) and ( |l4|) one can check that for the quark and 
squarks singlets the order of magnitude is given by (xD x ) ~ raw / Os- 
lo. Fig. (||) we show how B distributions change for different values of Mx- Each curve is the 
sum of the quark and squark singlets, gluon and gluino contributions. Note that d T D x > for all 
x because of the inhomogeneous term in Eq. (|?]). If the sparticle mass scale Msusy is around the 
electroweak scale, the final result depends feebly on Msusy [ 15 1 . 

The decay of X will produce neutral winos as well. The B and W° curves have similar shapes, 
the main difference between them being that, for a common scale Mx, W° FFs are always slightly 
larger than B FFs, since the coupling a*y is slightly stronger for W° than for B. The wino will 
eventually decay into the bino, the LSP. Its lifetime times its Lorentz factor is much smaller than 
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1 kpc/c, therefore neutral winos produced in the galactic halo or further will disintegrate into the 
LSP before reaching the Earth. The B spectrum on the Earth will be the sum of the B contribution 
produced at the decay spot of X plus the B contribution stemming from W° decay into B. 

The decay sequence of a neutral wino into a bino depends on the details of the scenario consid- 
ered. Even within the framework of the Constrained MSSM, the freedom is still large. Nevertheless, 
some benchmark points have been proposed [ 30 1 for study at the Tevatron collider, the LHC and 
e + e~ colliders. These points are consistent with different experimental constraints, as well as cos- 
mology, and can be regarded as generic in the whole parameter space. We shall be only concerned 



with the benchmark points A-D and G-M in Ref. |3C], where the LSP is mainly the bino and the 
next-to-LSP is the neutral wino, with some admixture of higgsinos. The two remaining points 
correspond to the so-called "focus-point" region at large scalar masses, that we shall not consider. 
In any case, our analysis covers a large region of the allowed CMSSM parameter space. 

The neutral wino decays mostly into sleptons and leptons, and it is followed by a decay of 
the slepton into the LSP. Therefore, only two-body decays are relevant. Given the generic case 
A — > B + C, we want to calculate which is the distribution of the decay product dnc/dx once the 
distribution for the decaying particle dnA_/dx is known. From the phase space of the disintegration 
process we obtain 

dn c {x) _ m\ pfmax dy dn A (y) 



d£ cr 2 Jy min y dy 



where we have defined 



Vmax = mm 



2m% X ' 



(20) 
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Umin — X ' (^) 

ax = m\ + rr? c — m 2 B , (22) 

cr 2 = ^X(m%m 2 B ,m 2 c ), (23) 

A(o, 6, c) = a 2 + 6 2 + c 2 - 2a6 - 2ac - 26c. (24) 
In addition there is the following kinematical constraint 

0<x <Jl^L (25) 
o"i - a 2 

The LSP flux reaching the atmosphere is the sum of the B contribution produced directly in 
the decay of X plus the B contribution produced by the decay of W° on its way to the Earth. The 
two-body decay of winos into binos pushes the momenta to lower values of x, making the direct 
B component dominate over the component produced by W° decay, for x > 0.1. 

So far we have ignored higgsino production, so some words are now in order. The magnitude of 
the gaugino and higgsino couplings to partons are comparable, hence their FFs will be also similar. 
For the benchmark points that we have analysed, higgsinos will decay on their way to the Earth 
into winos and subsequently into binos. Therefore, for the same kinematical reason as the wino 
case discussed above, the total number of binos on Earth coming from higgsino decay will be a 
small correction to the direct bino contribution, for x > 0.1. For the same reason, we have ignored 
binos coming from the decay of the charginos that are also produced in the parton shower. 

We compare in Fig. || the spectra of baryons, neutrinos and LSPs expected in the case that 
UHECRs are produced by the slow decay of a population of superheavy dark matter particles X 
clustered in the galactic halo. We take as dark matter particle mass Mx = 10 12 GeV and as scale 
at which SUSY switches on MgusY = 400 GeV. The baryon and neutrino curves are obtained from 
Ref. [15]. The neutralino spectrum shown corresponds to benchmark point C in Ref. [3C]. The final 



neutralino spectrum does not depend significantly on which benchmark point we choose, since the 
dominant contribution in the region of interest is given by the bino component produced directly 
byX. 



IV. CONCLUSIONS 

We have used the DGLAP formalism to calculate the spectrum of neutralinos produced by the 
decay of particles with mass Mx > 10 12 GeV. We have concentrated on the electroweak radiation 
of neutralinos by partons. This contribution has to be added to the on-shell contribution calculated 
with Montecarlo generators. We find that at large x the radiative contribution is slightly larger. 

For x > 0.3, LSPs dominate over baryons, photons and neutrinos and their flux on Earth 
is non- negligible. If UHECR are produced by the decay of superheavy particles and SUSY is a 
low energy symmetry of nature, then forthcoming observatories sensitive to weakly interacting 
UHECRs should detect a flux of ultra high energy neutralinos. 



ADDENDUM 



As this manuscript was being finished, Ref. [31] appeared. This paper also discusses neutralino 
production by the decay of super-heavy dark matter particles using DGLAP evolution. 
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FIG. 3. Neutrino (dashed lime), baryon (solid line) and neutralino (dotted line) spectra ex- 
pected from the decay of dark matter particles with Mx = 10 12 GeV clustered in the galactic 
halo. 
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